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Including the effect of thermal fluctuations in traditional computational fluid dynamics requires 
developing numerical techniques for solving the stochastic partial differential equations of fluctu- 
ating hydrodynamics. These Langevin equations possess a special fluctuation-dissipation structure 
that needs to be preserved by spatio-temporal discretizations in order for the computed solution to 
reproduce the correct long-time behavior. In particular, numerical solutions should approximate the 
Gibbs-Boltzmann equilibrium distribution, and ideally this will hold even for large time step sizes. 
We describe finite- volume spatial discretizations for the fluctuating Burgers and fluctuating incom- 
pressible Navier-Stokes equations that obey a discrete fluctuation-dissipation balance principle just 
like the continuum equations. We develop implicit-explicit predictor-corrector temporal integrators 
for the resulting stochastic method-of-lines discretization. These stochastic Runge-Kutta schemes 
treat diffusion implicitly and advection explicitly, are weakly second-order accurate for additive noise 
for small time steps, and give a good approximation to the equilibrium distribution even for very 
strong fluctuations. Numerical results demonstrate that a midpoint predictor-corrector scheme is 
very robust over a broad range of time step sizes. 



I. INTRODUCTION 



Modeling the effects of thermal fluctuations in spatially-extended systems requires introducing Langevin stochastic 
forcing terms in traditional deterministic models [TJI2]- Stochastic effects arise in fluid dynamics because of the random 
thermal motion of the molecules comprising the fluid at the microscopic level. Stochastic effects are important in 
flows at micro and nano scales typical of new nano- and micro-fluidic and microelectromechanical devices 0] , novel 
materials such as nanofluids [5], and biological systems such as lipid membranes, Brownian molecular motors, and 
nanopores [6]. Thermal fluctuations can be amplified by non-equilibrium effects and affect the macroscale, as in 
fluid mixing [5], propagation of fronts [SI [TO], combustion of lean flames, capillaries [ITJ[T2], and hydrodynamic 
instabilities [13-15 . Because they span the whole range of scales from the microscopic to the macroscopic [HH], 
fluctuations need to be consistently included in all levels of description, including continuum descriptions. 

Thermal fluctuations can be included in the classical Navier-Stokes-Fourier equations of fluid dynamics and related 
conservation laws through stochastic forcing terms, as first proposed by Landau and Lifshitz. The original formulation 
was for compressible single-component fluids jTBj. However, the methodology can be extended to other systems such 
as fluid mixtures [17] , chemically reactive systems [18] , magnetic materials [19] , and others [20] . The structure of the 
equations of fluctuating hydrodynamics can be, to some extent, justified on the basis of the Mori-Zwanzig formalism 
[211 122j . The basic idea is to add a stochastic flux corresponding to each dissipative (irreversible, diffusive) flux [2J, 
leading to a continuum Langevin model that ensures the correct equilibrium distribution. Specifically, statistical me- 
chanics tells us that the stationary (invariant) distribution at thermodynamic equilibrium is the Einstein distribution 
for isolated systems, and the Gibbs-Boltzmann distribution for systems in contact with a thermal bath. 

As model equations, here we consider the fluctuating Burgers equation in one dimension and the fluctuating Navier- 
Stokes equation in two and three dimensions. These stochastic conservation laws have non-dissipative (skew-adjoint) 
advective terms and dissipative (self-adjoint) viscous terms, as well as stochastic forcing terms that are in fluctuation- 
dissipation balance with the dissipative terms. Spatial discretizations of the corresponding continuum (functional) 
operators should preserve these (anti-)symmetry properties in order to obey a discrete fluctuation-dissipation balance. 
In general, constructing such spatial discretizations (coarse-grained Langevin models) is non-trivial and may be at 
odds with other considerations such as deterministic stability. For example, upwind discretizations commonly used for 
advection-diffusion equations add artificial dissipation to the equations, thus violating fluctuation-dissipation balance. 

Based on prior work by us and others, we construct spatial discretizations for the fluctuating Burgers and fluctuating 
incompressible Navier-Stokes equations that obey a discrete fluctuation-dissipation balance principle just like the 
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continuum equations. The main challenge is in constructing temporal integrators for the resulting large-scale system 
of stochastic differential equations. Ideally, the temporal integrators should have higher-order short-time accuracy 
but also lead to long-time dynamics that is in agreement with fluctuation-dissipation balance. Generalizing temporal 
integrators that are favored for Langevin equations in low-dimensional systems (e.g., molecular dynamics) [23] to the 
equations of fluctuating hydrodynamics would require implicitly handling the non-linear advective terms. Solving the 
resulting large-scale non-linear system of equations is computationally expensive in three dimensions for the Navier- 
Stokes equations. Mixed implicit-explicit Runge-Kutta schemes are used commonly in the deterministic method of 
lines, and can be extended to the stochastic setting with little effort, at least for the case of additive-noise equations 

EH- 

We derive the conditions for second-order weak accuracy of implicit-explicit predictor-corrector schemes for additive 
noise, and construct several candidate schemes. These semi-implicit schemes use an implicit midpoint rule for the 
diffusive and stochastic terms, which has the remarkable property that it gives the correct equilibrium distribution 
independent of the time step size. The non-linear advective terms are handled explicitly using an Euler predictor 
and a trapezoidal or midpoint corrector. We explain how to incorporate the stochastic forcing term in the resulting 
predictor-corrector schemes, and numerically study their performance on the model fluctuating Burgers equation. 
The numerical results suggest that the midpoint corrector is particularly robust. We then extend the spatio-temporal 
discretization to handle the incompressibility constraint present in the fluctuating Navier-Stokes equations in two and 
three spatial dimensions. We also include a passively- ad vected fluctuating scalar field, as encountered in practical 
problems. The schemes that we develop and analyze here have already been employed in Ref. |25j to construct a 
robust numerical solver for fluctuating incompressible flows, and to simulate the appearance of giant concentration 
fluctuations in diffusively mixing fluids [7] 18] . 

We begin by reviewing the fluctuating Navier-Stokes equations, and summarize a rather general formulation of 
finite-dimensional Langevin equations that obey a fluctuation-dissipation principle. In Section [IT] we explain in detail 
how to spatially discretize the fluctuating Burgers in one dimension so as to obtain a finite-dimensional system 
of Langevin equations with the proper structure. We turn our attention to temporal integrators for solving the 
resulting large-scale structured system of stochastic differential equations in Section |III| In Section |IIID| we design 
several implicit-explicit stochastic Runge-Kutta schemes that are second-order weakly accurate for additive noise and 
maintain fluctuation-dissipation balance even for large time steps. In Section |TV] we explain how the spatio-temporal 
discretization developed for the fluctuating Burgers equation can be generalized to the fluctuating Navier-Stokes 
equations with a passively-advected scalar field. The performance of the proposed schemes in the nonlinear (large 
fluctuation) setting is investigated in Section|Vj and some concluding remarks are given in Section VI Several technical 
calculations are detailed in the Appendix. 



A. Fluctuating Hydrodynamics 

The prototype stochastic partial differential equation (SPDE) of fluctuating hydrodynamics is the fluctuating Navier- 
Stokes equation. This equation approximates the dynamics of the velocity field v(r, t) of a simple Newtonian fluid in 
the isothermal and incompressible approximation, V • v = 0, 

p(d t v + v ■ V«) = -V7r + 77V 2 u + V • (2k B Tr])^ Z + / (1) 

where tt is the non-thermodynamic pressure, p is the (constant) density, ?/ = pv is the (constant) shear viscosity and 
v is the kinematic viscosity, and / (r,t) is an additional force density such as gravity [17] • Note that we prefer to use 
the standard physics notation instead of the differential notation more common in the mathematics literature, since 
the noise is additive and there is no difference between the different interpretation of stochastic integrals (e.g., Ito vs. 
Stratonovich) . In the momentum conservation law ([I]), the stochastic momentum flux is modeled using a white-noise 
random Gaussian tensor field Z (v,t), that is, a tensor field whose components are independent (space-time) white 
noise processes, 

{Z l3 {r,t)Z kl {r',t')) = (SikSj^Sit - t')S(r - r'). 

Note that in principle the stochastic momentum flux should have the symmetrized form (kBTrf) 1/2 (z + Z T ^j [T7]; 

however, for incompressible flow with constant viscosity this is not necessary |25j . 

The fluctuating Navier-Stokes equation, like other augmented Langevin equations of interest [26], obeys a 
fluctuation-dissipation principle, as explained more precisely in Section IB Specifically, (JlJ is constructed so that, 



at thermodynamic equilibrium, the invariant measure (equilibrium distribution) for the fluctuating velocities with 
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periodic boundaries is the Gibbs-Boltzmann distribution with a coarse-grained free energy or Hamiltonian given by 
the kinetic energy of the fluid, formally, 

P cq (v) = Z- 1 exp 

This is ensured by constructing the stochastic forcing term so that its covariance is proportional to the viscous 
dissipation operator 77V 2 . The advective operator (v ■ V), is, at least formally, Hamiltonian [27] in nature, which 
means that it preserves the equilibrium distribution dictated by the competition between the dissipative and the 
stochastic forcing terms. We argue that these well-known observations about the structure of the continuum equations 
should guide the construction of spatio-temporal discretizations [2"g] . 

We have formally written equation as an infinite dimensional stochastic differential equation. However, the 
interpretation of the nonlinear term v ■ V« requires giving a precise meaning to products of distributions, which 
cannot be defined in general and requires introducing some sort of regularization. An alternative is to define a 
discrete hydrodynamic field directly via some form of averaging of the molecular configuration of the fluid, and to 
obtain directly a finite-dimensional system of stochastic ordinary differential equations (SODEs) for the discrete 
variables through the Mori-Zwanzig formalism J55J I2H ISO]- While such an approach has certain advantages from 
a coarse-graining perspective, the notion of a continuum equation and the applicability of traditional methods for 
computational fluid dynamics is lost or at least obscured. 

Here we adopt a middle ground between the "continuum" and the "discrete" approach to fluctuating hydrodynamics. 
Specifically, we first spatially discretize the SPDE to obtain a system of SODEs, in the spirit of the "method of lines". 
Our focus here is on the temporal integrators for the resulting system of SODEs. We do not consider the convergence 
of the numerical method as the spatial discretization is refined, as one would in deterministic fluid dynamics. Rather, 
we fix the spatial discretization and assume that the hydrodynamic cells are sufficiently large, specifically, that they 
contain, on average, sufficiently many fluid molecules N p 3> 1. This ensures that the equilibrium fluctuations will be 

— 1/2 

on the order of 0(N P ) relative to macroscopic fields. We also assume that the transport coefficients have been 
renormalized to account for the finite number of fluid particles (molecules) used to define the hydrodynamic fields 
[H [3T] . There is strong numerical evidence that under these conditions spatio-temporal discretizations can correctly 
capture the leading-order (measurable) effects of fluctuations at large scales, such as fluctuation-driven transport in 
non-equilibrium systems [8], large-scale inhomogeneities arising during free fluid mixing [25], and diffusive effects on 
the very long-time dynamics such as drifts in propagating fronts [35] and shocks [33] , 
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B. Fluctuation-Dissipation Balance in Generic Langevin Equations 

The fluctuating hydrodynamic formalism finds its foundation in the theory of coarse-graining [34 . One of the 
central objects in the theory is the coarse-grained Hamiltonian or a coarse-grained free energy, which determines the 
Gibbs-Boltzmann equilibrium probability density for the coarse variables x € K. , 



P eq (x) = Z 1 exp 



H(x) 



(2) 



where Z is a normalization factor and the scaled temperature k B T sets the unit of energy. We will set fcgT = f 
for simplicity through the remainder of this paper. Note that in many cases the actual Hamiltonian (total energy) 
may be expressible in terms of the coarse-grained variables and is strictly conserved [5], but this is not the case for 
isothermal systems. 

As discussed at length by Grabert [35] . a Markovian approximation within the Mori-Zwanzig formalism can be used 
to obtain a generic Langevin equation of motion |26j for an arbitrary choice of the microscopic ensemble, 

d t x = -N (x) °— + (2k B T) 1/2 B (x) W(t) + (k B T) |- • N* (x) , (3) 

where W € M. Nw denotes white noise, the formal temporal derivative of a collection of independent Brownian motions, 
and an fto interpretation is assumed. We will typically suppress the explicit dependence on x and write the mobility 
operator as N = N (x). In the generic Langevin equation ([3]), the skew-adjoint operator 

S = -S* = - (N* - N) 



4 



generates the "conservative" part of the dynamics, and the self-adjoint positive semi-definite operator 

M = M* = - (N + N*) t 



generates the "dissipative" part of the dynamics, since 



dH 
~dt 



dHy BH 

~E — N— = -Re 

Ox I ox 



dH_\ 
~dx) 



M 



dH 
dx 



< 0. 



The operator B (x) is constrained, but not uniquely-determined, by the fluctuation- dissipation balance condition 

BB* = M. 

The last term in (|3| is an additional "spurious" or "thermal" drift term whose form depends on the particular interpre- 
tation of the stochastic equations [31)], which we take to be in the sense of Ito. In the case of the fluctuating Burgers 
and Navier-Stokes equations this term will vanish. 

The dynamics (|3| is ergodic and time reversible (under an appropriate parity transformation for the variables) with 
respect to the distribution B. It is not hard to show using the corresponding Fokker-Planck Equation that ([2| is an 
equilibrium distribution for O), as desired, at least if one makes the nonrestrictive assumption that M >~ 0. More 
precisely, the invariant measure for the coarse-grained variables is dfJ, eq = /i cq (dx) = P cq (x) dx. We will assume here 
that the dynamics is ergodic with respect to a unique equilibrium distribution, which may not be the case if H (x) is 
not defined everywhere. It is important to point out that when Langevin equations of the form (|3| are applied in a 
non-equilibrium setting, the dynamics may not be ergodic with respect to any distribution, even at steady state. 



II. FLUCTUATING BURGERS EQUATION 

The analysis and numerical solution of the incompressible Navier-Stokes equation is complicated by the presence 
of the incompressibility constraint. We begin our discussion by constructing a spatial discretization of the simpler 
unconstrained fluctuating Burgers equation for the random field u(x,t), 

9 - 

d t u + cu d x u = v d xx u + (2v) 2 d x Z, (4) 

where v is a diffusion coefficient and c sets the scale for the advection speed. This equation mimics some of the 
properties of the fluctuating Navier-Stokes equation ([I]), in particular, it obeys a fluctuation-dissipation balance 
principle with respect to the Gibbs-Boltzmann distribution with a Hamiltonian H — j dxu 2 /2. The fluctuating 
Burgers equation can also be written in conservative form 



d t u = -cL 



c— — vd x u — (2v) 2 Z 



showing that the total momentum J u dx is conserved with periodic boundary conditions. Note that here the stochastic 
forcing term is linear and involves the spatial derivative of white noise [37]: rather than white noise itself as in the 
stochastic Burgers equation studied, for example, in Ref. [38 . Equations of this type arise as coarse-grained models 
of the behavior of one dimensional lattice gases, such as the asymmetric excluded random walk model [37 . 

In this section we show how the fluctuating Burgers equation can be spatially discretizcd in a manner that leads to 
a generic Langevin equation of the form (J3j) . This construction will be extended to the Navier-Stokes equations with a 
passively-advected scalar in Section |IV| Our approach to the spatial discretization follows standard practice in deter- 
ministic fluid dynamics. Specifically, we construct the spatially-discrete system by combining locally-accurate spatial 
discretizations of the differential operators (e.g., gradient, divergence and Laplacian) that appear in the the SPDE. 
However, in addition to focusing on accuracy and stability when choosing the spatial discretization, we pay particular 
attention to preserving fluctuation- dissipation balance. This means that we want to obtain a system of SODEs whose 
structure is given in ^ and whose invariant distribution (equilibrium distribution) is a natural discretization of the 
Gibbs distribution dictated by equilibrium statistical mechanics. 



A. Continuum Fluctuating Burgers Equation 

One can, at least formally, consider a generic Langevin equation for an infinite dimensional field [2 . The fluctuating 
Burgers equation Q is a prototype of such an equation. In this formalism the coarse-grained Hamiltonian is a 
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functional of the field and the partial derivatives should be interpreted as functional derivatives, and contractions by 
a field imply integrations over the spatial domain. For equation Q , the (formal) free energy functional is 

r u 2 

H[u(x,t)} = I -dx, (5) 

so that 

dH _ SH[u(x,t)\ _ 
du Su 

The dissipative and fluctuating dynamics in Q are generated by the constant operators, 

M = -vd 2 xx and B = v^d x , 

which in higher dimensions become multiples of the Laplacian and divergence operators, respectively. The conservative 
dynamics for the Burgers equation is Hamiltonian and generated by the skew-adjoint linear operator S (u) defined 
through its action on a field w (x, t) [27\, 

S (u) w = — - [ud x w + d x (uw)] . (6) 
o 

The v ■ Vi> term in the higher-dimensional fluctuating Navier-Stokes equation ([I]) can similarly be written in terms 
of a skew-adjoint operator, although there are some complications in handling the divergence-free constraint |39) . 

A detailed description of the meaning and importance of the Hamiltonian nature of the nonlinear deterministic 
dynamics and the Poisson bracket associated with S is given in Refs. [HHTJ. F° r our purposes, the most important 
property of Hamiltonian dynamics is that it is incompressible in phase space, 

». SW _|-.jV(.)-o. (7) 

This implies that the dynamics of the inviscid Burgers equation preserves not just functions (such as the Hamiltonian 
itself) but also phase-space measures (such as the Gibbs distribution), and thus any probability density that is 
a function of H only is a candidate equilibrium distribution. The inviscid Burgers equation may also be written 
in Hamiltonian form using the Hamiltonian H = J (u 3 /6) dx with S = —d x [27J HO]- However, in fluctuating 
hydrodynamics the choice of the coarse-grained Hamiltonian is dictated by statistical mechanics and the equilibrium 
Gibbs distribution is maintained via the fluctuation-dissipation balance between the viscous and stochastic terms. 



B. Discrete Fluctuating Burgers Equation 

The preceding discussion of how the fluctuating Burgers equation can be written in the form of a generic Langevin 
equation ^ is formal and merely informs our choice of spatial discretization. The discretized u — {m, . . . , u^} can 
be thought of as a finite-volume representation of the field u{x,t) on a regular grid with spacing Ax, specifically, Uj 
can be thought of as representing the average value of u(x, t) over the interval (cell) [jAx, (j + I) Ax]. As we already 
discussed, this is merely a formal association and the actual physical object is the discrete (coarse-grained) u (t) and 
not the hypothetical u(x,t). Similarly, the spatially-discretized collection of white noise processes can 
formally be associated with the space-time white noise Z. 

We take the coarse-grained Hamiltonian function to be the natural (local equilibrium [53]) discretization of ([5]), 

ff(«)=Ev«i. (8) 

We will construct a spatial discretization that leads to a finite-dimensional generic Langevin equation of the form ^ 

e,u = s d £ + JLD™ + (Z) 1,2 D lW{t) . (9) 



du Ax du \Ax 

Here W is a vector of N w independent white-noise processes (formally, time derivatives of independent Wiener 
processes), TD\ is a matrix representing the spatial discretization of the divergence operator, such that L>2 = D\D\ 



() 



is a symmetric negative-semidefinite discretization of the Laplacian operator. This system of SODEs has as an 
invariant distribution the Gibbs distribution ^ if S is an antisymmetric matrix discretizing (|6| that satisfies 



d_ 
du 



S{u) 



E 



dS 



for all k. 



(10) 



We now construct specific finite-difference operators for D\ and S. 

A particularly simple choice that also generalizes to higher dimensions [41] is to associate fluxes with the half-grid 
points (faces of the grid in higher dimensions), and to define 



(£>iW) ? 



giving {D\u) 



Uj+l 



This construction gives the familiar three-point discrete Laplacian (2c? + 1 points in dimension d), 

Uj—i — 2uj + uj + i 



Ax 2 



(11) 



and is therefore an attractive choice that satisfies the discrete fluctuation-dissipation principle [41] . If periodic bound- 
ary conditions are imposed, we set Uq = itjv and mjv+i = i*i and Wi = W N+ i (i.e., iV^, = N). For Dirichlct boundary 
conditions we fix uo and ujv+i at specified values and do not need to impose any boundary conditions on W (i.e., 
N W =N + 1). 

— — * 

A natural choice for S is formed by choosing a skew-adjoint discretization D\ = —D 1 of d x , in general different 
from Di, and discretizing ^ directly as 



(Su) 1 = -- 



(S 1 u) j + (Bxti 9 ) 



where m 2 = 



, ttfy}. We choose £>i to be the second-order centered difference operator 



D lU ) = 



2Aa 



leading to an explicit expression that makes it clear that Su is a discretization of —cuu x , 



2Ax 



2 2 



2Ax 



+ Uj + Uj + i ^ ^ Uj+i — Uj-i 
2Ax 



The above discretization of the advective term has been considered frequently in the literature, as discussed in detail 
in Ref. [10]. It can be seen as a weighted combination of the "convective" and the "conservative" forms of advection 
[42] with weights 1/3 and 2/3, which is the unique choice of weights that gives a conservative and skew-adjoint 
discretization of advection. It is important to note that one can write the nonlinear term in conservative form, 



(Su). = - t 



'it*, -u 2 
Ax 



wher" ' 



u 2 + UjUj + i + u 2 +1 



J+5 



(12) 



Due to the skew-symmetry, in the absence of viscosity the total "energy" Q is conserved for periodic systems. It can 
also easily be shown that that the condition ([7]) is satisfied and therefore this particular discretization of the advective 
term preserves the Hamiltonian structure of the equations |40| . 

Putting the pieces together we can write the semi-discrete fluctuating Burgers equation as a system of SODEs, 
j = l,...,N, 



duj 
~dt 



QAx 
v 



(uj-i + Uj + Uj+i) (uj+i - Uj-i) 



(13) 



+ £2 K'-i - 2«i + uj+i) + (t) - W i4 (t)) 



With periodic boundary conditions, this stochastic method of lines 
energy |8]) and the total momentum 



discretization strictly conserves the total 



i(u) = Ax Uj 
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The equilibrium distribution is the discrete Gibbs-Boltzmann distribution 

N 



Ax ^ 2 



P cq («) = Z- 1 exp 

where mo is the initial value for the total momentum. In the next section, we construct efficient temporal discretizations 



2 ^ * 



5 Ax Uj — m 



of ( 13 ) that preserve these properties as well as possible. 



III. WEAKLY-ACCURATE TEMPORAL INTEGRATORS 

In this section, we consider weak temporal integrators for coupled systems of stochastic ordinary differential equa- 
tions (SODEs) that arise after spatial discretization of the equations of fluctuating hydrodynamics, such as the system 



(13). To make the discussion more general and applicable to a large range of generic Langevin equations, we consider 

! differential equations 

a (x) + KW (i) = [L (x)] x + g (x) + KW (t) . (14) 



the system of nonlinear additive-noise differential equations 

dx 
~dt 



There are N v independent variables x (t), and W (t) denotes a collection of N w independent white-noise processes, 
formally identified with the time derivative of a collection of independent Brownian motions (Wiener processes) . Here 
If is a constant matrix, and we have used the more natural differential notation since there is no difference between 
the different stochastic interpretations (e.g., Ito and Stratonovich). For fluctuating hydrodynamics applications, we 
split the drift a (x) into a diffusive term [L (x)] x and an advective term g (x). This becomes particularly important 
when considering semi-implicit temporal discretizations since diffusion often needs to be treated implicitly for stability 
reasons. In general, L(x) may depend on x since the transport coefficients (e.g., viscosity) may depend on certain 
state variables (e.g., density). 

Multiplicative noise poses well-known difficulties with constructing higher-order temporal integrators |43j . There 
are many important fluctuating hydrodynamic equations with additive noise, such as the fluctuating Navier-Stokes 
equation ([!]), in which L (x) = L = const is a multiple of the (discrete) Laplacian and K (x) = K = const, is a 
multiple of the (discrete) divergence. For generality in Appendix [A] we also consider constructing first-order accurate 
schemes for the case of multiplicative noise, where K (x) also depends on x, and here we focus on second-order 
accuracy for additive noise. 

There is an extensive literature on numerical methods for finite-dimensional systems of SDEs. At the same time, 
efficient numerical solution of the types of systems of SDEs appearing in the stochastic method of lines [H] for fluctu- 
ating hydrodynamic equations limits the choices of practicable techniques. One of the most important characteristics 



of systems such as (13) is the presence of a large number of length scales and associated relaxation times (hydro- 
dynamic modes). This intrinsic stiffness is particularly prominent for diffusive terms and makes the construction of 
stochastic integrators particularly challenging. A powerful class of integrators for systems with a broad spectrum of 
relevant time scales are exponential integrators [44] [45] . These integrators apply a local linearization [46] and use 
the matrix exponential exp (LAt) to capture the dynamics at a time step At that under-resolves the stiffest modes. 
Unfortunately, the computation of the matrix exponential is only practicable when it is simple to diagonalize L, as is 
the case when a Fourier basis is used for periodic systems |47j . Multi-step schemes are widely used in the deterministic 
context because of their favorable stability properties and reduced computational effort per time step [48] [49] . Weak 
multistep schemes have not received much attention in the literature, and their analysis requires developing novel 
techniques that are outside the scope of this work. 

With these considerations in mind, we focus here on one-step Runge-Kutta schemes. They have the advantage 
of potentially attaining high order of accuracy without requiring evaluation of derivatives of a (x). We will consider 
numerical schemes that only require the generation of N w = O (N v ) Gaussian random variables and do not require the 
solution of nonlinear systems of equations. While this limits the robustness of the schemes in the nonlinear setting, 
it is the only type of method that is practicable for the sort of large systems of SODEs that arise when discretizing 
hydrodynamic SPDEs. First we discuss the simpler case of explicit stochastic Runge-Kutta integrators, and then we 
consider stochastic variants of two-stage implicit-explicit Runge-Kutta (IMEX-RK) [50] integrators. 

A. Conditions for Second-Order Weak Accuracy 

In this section we discuss second-order weakly accurate one-step temporal discretizations with a constant time step 



At, and denote the numerical approximation x n w x (nAt) . The most fundamental temporal integrator for (14) is 



the weakly first-order accurate Euler-Maruyama scheme, 

x n+1 = x n + Ata n + At^K n W n , (15) 

where the superscript denotes the time level at which the term is evaluated, for example, a" = a(x n ). Here W n is 
a vector of m independent standard Gaussian variates (i.e., normally-distributed pseudorandom numbers with mean 



zero and unit variance), generated independently at each time step. The SDE ( 14 ) can, in fact, be defined through the 



limit At — > of the scheme (15). The stochastic term At 2 W n represents the (Wiener) increment of the underlying 



Brownian motion s oy er the time step. Alternatively, one can view (15) as an application of the deterministic explicit 



Euler method to (14), with the discrete white noise At 1 / 2 W n representing the rough forcing W (t). This viewpoint 



is particularly useful when extending higher-order standard deterministic schemes to the stochastic context. 



Let us consider one-step schemes for the general nonlinear additive-noise system of SDEs ( 14 ) . The general theory 
of weak accuracy for stochastic integrators is well-established and reviewed, for example, in Section 2.2 of [5T]. The 
key result is that, under certain assumptions, second-order weak accuracy is achieved if the first 2-2 + 1 = 5 moments 
of the numerical increment Ax n = x n+1 — x n match the moments of the increment x (nAt + At) — x (nAt) to order 
O (At ) . The required moments can be obtained from the well-known weak expansion 



1 / \ At 

x a (nAt+At) = x a {nAt)+At?K af} W$+Atal + - [At 2 a^ + At^K^W?) (d 7 a^) + — K lf K Se (d 7 d s a^) + 0(At 5 / 2 ), 

(16) 

where a repeated index implies summation and the short-hand notation d~ ( = d/dx^ is employed. This expansion is not 
directly useful for numerical approximations since it requires evaluating derivatives of the drift function. Instead, we 
employ Runge-Kutta schemes and then ensure second-order weak accuracy by matching the moments of the numerical 
increments to (16). The details of these calculations are given in Appendix [A| 



A well-known weakly second-order accurate predictor-corrector scheme that is consistent with the conditions derived 
in Appendix [A] is a two-stage explicit trapezoidal method [351 H3J [52] ■ In this scheme the first stage is an Euler- 
Maruyama predictor, and the corrector stage is an explicit trapezoidal rule, 



x n+1 = x n + Ata n + At?KW 7 \ 

x n+1 = x" + ^(a n+1 +a n ) + At^KW n . (17) 

One can naively obtain this scheme from the classical deterministic predictor-corrector algorithm by thinking of 
F n = At~ x l 2 K.W l as a constant applied forcing. The fully explicit nature of this method severely restricts the 
time step due to stability limits arising from the stiff diffusive terms in fluctuating hydrodynamics. We will consider 
semi- implicit RK schemes in more detail in Section [HI D| 



B. Equilibrium Fluctuation Spectrum 



An important property of Langevin-type equations, including those of fluctuating hydrodynamics, is the existence of 
a non-trivial stationary distribution (invariant measure) . It is important for numerical schemes to have an equilibrium 
distribution that is in, some appropriate sense, close to that of the stochastic differential equations. A recently 
proposed-approach 53J is to add a Metropolis-Hastings acceptance-rejection rule to a classical integrator such as 
the Euler-Maruyama scheme. This "Metropolization" ensures that the equilibrium distribution of the numerical 
approximation is controlled; however, this is done at the cost of reducing the temporal accuracy because of rejections. 
It is therefore important to ensure that the non-Metropolized numerical scheme produces a good approximation to 
the equilibrium distribution, so that rejections are infrequent. 

Mattingly et al. |54] show that in some appropriate metric the invariant measure (which is assumed to exist) of the 
numerical scheme has the same order of accuracy as the weak order of accuracy over finite time intervals. This only 
provides an asymptotic error bound, however, and does not provide an estimate of the actual error. By focusing on the 
linearized equations of fluctuating hydrodynamics one can easily obtain explicit estimates for the invariant measure 
of a given numerical scheme and thus understand the nature of discretization errors in the long-time dynamics. This 
approach was used by some of us in an earlier publication |41j to analyze and improve explicit Runge-Kutta schemes 
for compressible fluctuating hydrodynamics. Here we briefly review the main results and discuss some generalizations. 

We consider the linear system of additive-noise SDEs 



— = Lx + KW (t) 



(18) 
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A general linear one-step temporal scheme for this equation has the form 

x n+1 =Qx n + AtiRW n , 

where Q and R are some iteration matrices. Since this is a linear equation forced by a Gaussian process, the solution 
is a Gaussian process. The equilibrium or steady-state covariance C \ t = (x 11 ( xn )*) of this linear iteration is the 
solution of the linear system (see, for example, the derivation in [41] ) 

QC At Q* - C At = -At RR*. (19) 

In the limit At — > any consistent and stable numerical scheme should give the correct equilibrium covariance 
C = (x (t) x* (t)), which is the solution to [Tl |4"TI 155] 

LC + CL* = KK*. (20) 



Equation ( 18 \ can easily be solved explicitly to obtain an exact exponential integrator for which Q = exp (LAt) 
follows from the deterministic variation-of-constants formula. This exponential scheme will be an exact integrator for 
([18} if and only if 

RR* = A^ 1 [C - QCQ*} = At- 1 [C - exp (LAt) C exp (L* At)] . (21) 

In general, one cannot write an explicit solution to this equation unless one can explicitly diagonalize L and C in 
some basis. 



Implicit Midpoint Rule 



Runge-Kutta schemes approximate the matrix exponential exp (LAt) with a polynomial (for fully explicit schemes) 
or a rational (for semi-implicit schemes) approximation. An important example is provided by the implicit midpoint 
(equivalently, trapezoidal) method (Crank- Nicolson scheme) applied to the linear problem (181, 



x n+1 = x n 



At 



L{x n + x n+1 )+At 1 2KW n . 



In this scheme the iteration matrix Q is a 1 — 1 Pade approximation of the matrix exponential, 



Q 



LAty 1 



LAt\ 



= exp (LAt) + O (At 3 ) 



(22) 



(23) 



and R = (I — LAt/2) 1 K. It is not hard to show that the implicit midpoint scheme leads to the correct equilibrium 
covariance C for any time step size since 

RR* = At 1 [C - QCQ*} , 

as seen from a straightforward explicit calculation, 
At" 1 [C-QCQ*] 



Ar 1 I 



LAt 



-1 r 



LAt\ 



CI 



L*At 



LAt 



CI 



L*At 



LAty 1 



(-LC-CV) [I- 



L*At 



KK* 



2 

L* At 
2 



L* At 

2 

RR*. 



An alternative derivation of the fact that ( |22| ) gives the correct steady-state covariance for any time step size At can be 
found in the Appendix of Ref. [25] . That derivation is based on showing that the iteration ( 22 ) is a Metropolis-Hastings 
Monte Carlo algorithm to sample the invariant distribution of (18). 

The implicit midpoint rule can easily be generalized to the nonlinear system (14), 

„n+l I „ 



At a 



AtiKW n , 
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which can be shown to be weakly second-order accurate. This scheme is a particularly good candidate for solving 
Langevin-type equations because it is a time-reversible and quasi-symplectic integrator [23] that exactly conserves 
all quadratic invariants (e.g., a quadratic Hamiltonian). However, it requires the solution of a nonlinear system of 
equations at every time step. This nonlinear system of equations may not have a unique solution and is in general 
too expensive to solve for large-scale hydrodynamic calculations. In the special case of the stochastic Burgers or 
Navier-Stokes equations, the only nonlinearity in a (x) comes from the advective term, which has the special form 

(to within irrelevant constants) [S (x)] x and can be linearized as S (x n+ ^ x, where x n+ ^ is a mid-point estimate 

that has to be obtained via a predictor stage. Such an approach gives a scheme that only requires solving a linear 
systems in each time step, while still preserving quadratic invariants (e.g., total kinetic energy). It is, however, not a 
time-reversible scheme. Furthermore, solving the non-symmetric systems that arise when advection is discretized in a 
semi-implicit manner poses a significant linear algebra challenge, especially when constraints such as incompressibility 
are included. For this reason, in the next section we consider implicit-explicit Runge-Kutta schemes in which only 
diffusive terms are, potentially, discretized implicitly. 



C. Fully Explicit Runge-Kutta Scheme 

We now illustrate how the conditions derived in Appendix [A] can be used in practice to construct a fully-explicit 
three-stage Runge-Kutta (RK3) integrator. In Ref. [33J , an algorithm for the solution of the compressible equations 
of fluctuating hydrodynamics was proposed, which is based on a well-known three-stage, low-storage total variation 
diminishing Runge-Kutta (RK3) scheme 56J . The RK3 scheme is a simple discretization for the deterministic com- 
pressible Navier-Stokes equations that is stable in the inviscid limit, even when slope- limitcrs are omitted for the 
convective terms. A stochastic version of the RK3 scheme was analyzed in Ref. [JT] based on a linearized analysis, 
and employed in Ref. [35] along with a staggered spatial discretization. In Ref. [H], nonlinearities were not considered 
and therefore the scheme presented there is not second-order weakly accurate for nonlinear equations. 

Each time step of the RK3 algorithm is composed of three stages, the first one estimating x at time t = (n + I) At, 
the second at t = (n + \)At, and the final stage obtaining a third-order accurate estimate at t = (n + l)At. Each 
stage consists of an Euler-Maryama step followed by weighted averaging with the value from the previous stage, 



x n+1 



--x n + Ata n + dd?K{axW n A + PiW n B ) 



x 



a 3 . 

2 — T» 

x n+1 =\x r 



x n+1 + At a n+1 + At*K{a 2 W n A + (3 2 W n B ) 
x n+ ^ +Ata n+ ^+A^-K(a 3 W 1 A :+f3 3 W B ) 



(24) 



Here ~W\ and W B are two independent vectors of i.i.d. normal random variates that are generated independently 
at each RK3 step, and the weights a and f3 are to be determined. In principle one could use a third sample Wq; 
however, this increases the cost of the method and is insufficient to yield a weakly third-order accurate scheme. 
Following Ref. [ST], we can fix a 3 and j3 3 by making the arbitrary choice that after all the stages are combined the 
stochastic increment in x n+ 1 be AtiKW n a . 

Runge-Kutta schemes of the above form have been analyzed in Ref. |57] and the moment conditions for weak 
accuracy of any order derived. Up to second-order one can easily obtain these conditions by explicit Taylor series 
expansion and comparison of the moments of the numerical increment to those in (16 1. This gives two quadratic 
equations for the weights a and j3. Two more equations can be obtained by asking for third-order accuracy of the 
static covariance C^t — C + At 3 AC + 0(At 4 ) in the linear case. This simple calculation consists of applying the 
RK3 scheme to the linear equation ( |18[) in order to extract an explicit expression for Q and R, and then substituting 
these expressions in (19) and using (|2*0"|) to eliminate KK*. By equating the coefficients in front of terms of lower 
order in At, we can ensure that the error in the stationary covariance is of order At 3 . In particular, equating the 
coefficients in front of the terms involving L 3 C and L 2 CL* to zero gives two additional quadratic equations for the 
weights a and /3. The solution of the resulting system of four equations for the weights oti, tijft and /3 2 gives a 
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stochastic RK3 scheme that is of weak order two in the general case and of order three in the linear case, 

ot\ =012 = (X3 = 1 
(2 V2± y/3) 



5 

„ (-4V2±3V3) 
ft- g 

_[72t2%/3] 
k - 10 ' 

For fluctuating hydrodynamics applications, we recommend the upper sign since it gives better discrete static structure 
factors for a model one-dimensional stochastic advection-diffusion equation (see Ref. [H] for an illustration of this 
type of calculations) . 



The RK3 scheme (24) suffers from a severe time step limitation due to the explicit handling of the diffusive terms. 



We consider alternative semi-implicit methods next. 

D. Semi-Implicit Runge-Kutta Temporal Integrators 

In this section, we construct two-stage second-order implicit-explicit Runge-Kutta schemes for solving a system of 
SDEs 

dx 

— = L(x)x + g(x) + KW (i), 
at 

where g (x) denotes all of the terms handled explicitly (e.g., advection or external forcing). We will consider schemes 
that at time step n require solving only linear systems involving the matrix L n = L(x n ), and are parametrized by 
a vector of weights w (Butcher tableau). The first stage in any of these schemes is a predictor step to estimate 
x rj x (nAt + lOaAt), where w 2 is some chosen weight (e.g., w 2 = 1/2 for a midpoint predictor). The corrector 
completes the step by estimating x n+1 at time (n + l)At, 

x = x n + (w 2 - Wl )AtL n x n + Wl AtL n x + w 2 Atg n + (w 2 At)^ KW™, 

x n+l = x n + (1 - w 3 - w 4 )At L n x n + w 3 At L n x + w 4 At L n x n+1 

+ w 5 At(L-L n )x + w 5 Atg + (l-w 5 )Atg n 

+ (w 2 At^ KW'l + {{I - w 2 ) AtY KW%. (25) 

1/2 

The intuition behind the handling of the stochastic increments in the predictor/corrector stages is that (w 2 Ai) 1/z Wl 

1/2 

samples the increment of the underlying Wiener processes over the time interval w 2 At 1 while ((1 — w 2 ) At) W 2 
represents the independent increment over the remainder of the time step. 

According to the calculations detailed in Appendix [S] in order to be second-order weakly accurate in the case of 
additive noise, the weights w should satisfy the conditions 

w 2 w 5 = ^, w 2 w 3 + w A = ^. (26) 

We now catalog some possible choices of these weights that satisfy these conditions. We will test the performance of 
these schemes on the fluctuating Burgers and Navier-Stokes equations in Section [Vj 

1. Explicit Midpoint Predictor- Corrector Scheme 

If wi — and w 4 — we obtain fully explicit schemes. In this case w% = w§ and therefore the splitting into implicit 
and explicit parts does not matter, and the only function that appears in the scheme is a(x) = L(x)x + g(x). The 
only parameter of choice is w 2 . A well-known method that fits the above format is th e ex plicit trapezoidal predictor- 



corrector method, w 2 — 1, and therefore w$ = w 3 = \. This is exactly the scheme (17 1. For linear equations, this 
scheme gives static covariances accurate to second-order, Ca* = C + At 2 AC + 0(At d ). A notable advantage of this 
scheme is that it requires generating only a single random increment per time step. 
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In order to choose the "best" u>2, we can look at the accuracy of the static covariance in the linear case, just as 

A simple calculation shows that in order to obtain third-order accuracy 
0(At 4 ), we need to use a mid-point predictor stage, w 2 = 1/2 and 



we did for the RK3 scheme in Section III C 
of the static covariance, C/^t = C + ApAC 
therefore w$ = W3 = 1. This gives the explicit midpoint predictor-corrector method 



„"+l 



At 



At a 



At 
~2 



At 



(27) 



which requires generating two random increments per time step. It is important to point out that for advection- 
diffusion problems the time step At needs to be substantially smaller than the stability limit for the explicit midpoint 



scheme (27) to realize its asymptotic order of accuracy and give substantially more accurate static covariances than 



the explicit trapezoidal scheme (17) 



2. Implicit Trapezoidal Predictor- Corrector Scheme 



We now extend the implicit midpoint scheme (22) to the case when some terms, such as diffusion, are handled 
implicitly. In Ref. |25j the fluctuating Navier-Stokes equation was solved using a scheme for which w\ — 1V4 — w$ — 
1/2, W2 = 1, W3 = 0, giving the implicit trapezoidal predictor-corrector method, 

x n+1 = x n + ^L n (x n + x n+1 ) + Atg n + At^KW n 

x n+i = ^« + ^ i »( K « + ^+i) + ^(Z" +1 _ Iy ») 5; »+i + ^(g»+i +9 n) + Aa J ft:vr". (28) 
If g = this scheme is stable for any time step, more precisely, it is A-stable (see Appendix [b] for L-stable schemes) . 



The scheme has the great advantage that for the linearized equation ( 18 ) it gives the correct stationary covariance 
regardless of the time step, while only requiring a single random increment per time step. We emphasize that the 
dynamics of the fluctuations is not correctly reproduced for large time step sizes, as we discuss further in Appendix 

m 



3. Implicit Midpoint Predictor-Corrector Scheme 



In the scheme (28), a trapezoidal approximation is used for the explicit fluxes, just as in (17). Another candidate is 



to use a midpoint approximation for the explicit fluxes, as in (27), obtained by using u>2 = 1/2 and therefore W5 = 1. 
We also choose a corrector stage where the implicit part is the implicit midpoint rule (22), which requires choosing 
W3 = and u>4 = 1/2. We are left with a choice for W\ in the predictor stage, 



:c' M " : : x" (\~wi) L'V'A/ ! l r l L"i:" + s - 



„«+l 



At 



At 



-L n (x n + x n+1 ) 



At(L ' 



L n )x n+ i + Atg n+ ; 



At 



(29) 



Two obvious choices are w\ — 1/2, which for linear equations makes the predictor stage a backward Euler step with 
time step size At/2. An alternative is to use W\ = 1/4, which for linear equations makes the predictor stage an 
implicit midpoint step with time step size At/2. For the linearized equation ( |18[ ) the predictor stage does not actually 
matter since it is only used in evaluating the nonlinear terms in the corrector stage. Therefore, we will compare the 
implicit midpoint schemes with w\ = 1/4 and w\ = 1/2 numerically in Section |v| 



IV. FLUCTUATING NAVIER-STOKES EQUATION 

The implicit-explicit schemes discussed in Section |III D| are general schemes suitable for unconstrained SDEs and 
cannot directly be applied to the fluctuating Navier-Stokes equation (fl]). Some care is required in handling the 
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incompressibility constraint in a computationally-efficient manner without compromising the stochastic accuracy. 
The spatio-temporal discretization we analyze here was proposed and applied in Ref. |25j ; here we provide additional 
analysis and a discussion of alternative approaches. 



A. Continuum Equations 



In principle, the incompressibility constraint can be most easily handled by using a projection operator formalism 
to eliminate pressure from ([!]) and write the fluctuating Navier-Stokes equation in the form 



d t v = V 



Vv + vV 2 v + (2iyp- 1 k B T) 2 V ■ Z t 



(30) 



Here "P is the orthogonal projection onto the space of divergence-free velocity fields, "P = I— Q {'DQ) 1 X> in real space, 
where D = V- denotes the divergence operator and Q = V the gradient operator with the appropriate boundary 
conditions taken into account. We only consider periodic, no-slip and free-slip boundaries. With periodic boundaries 
we can express all operators in Fourier space and "P = I — k~ 2 (kk*), where k is the wavenumber. The application 
of the projection to the right hand side ensures that V • v = at all times if the initial condition is divergence free. 
The divergence-free constraint is a constant linear constraint and the projection restricts the velocity dynamics to 
the constant linear subspace of divergence-free vector fields. The projection operator can be applied in more general 
settings, notably, in cases where the constraints are nonlinear and the noise is multiplicative; however, the resulting 
expressions are rather complex especially in the stochastic setting [351 US] ■ 

In practice, the fluctuating velocities modeled by (30) advect other quantities, and it is this coupling between 
the velocity and other equations that is of most interest. Perhaps the simplest example is provided by a stochastic 
advection-diffusion for the concentration or density c{r,t) of a large collection of non- interacting passive tracers. 
For example, c(r,t) might corresponds to the light intensity pattern of ffuorescently-labeled molecules suspended 
in the fluid in a Fluorescence Recovery After Photobleaching (FRAP) experiment. In general, the equation for the 
concentration has multiplicative noise [25]. This arises because the coarse-grained free energy functional H [x (r, t)] = 
H(v, c) includes a contribution from the entropy of the passive tracer which is, in general, a non-quadratic even 
if local functional of c. For illustration purposes we can take a separable quadratic Hamiltonian (i.e., independent 
Gaussian fluctuations in velocity and concentration), 



H (v, c) = H v (v) + H c (c) 



and write the the model additive-noise tracer equation 



v dr 



k B T 
2e 



z 2 dr, 



d t c ■ 



-v ■ Vc + xv c + V 



(2ex) 5 ^ c 



(31) 



The multiplicative noise case is not considered herein. Physically, e measures the degree of coarse graining, e ~ Np 1 , 
where N p is the number of tracer particles per coarse degree of freedom. Note that (31 ) is a conservation law because 
v ■ Vc = V • (cv) due to incompressibility. 



The coupled velocity-concentration system ( 30|31 ) can formally be written in the form The chemical potential 
fi (c) = dH/dc ~ c. The mobility operator can be written as a sum of a skew-adjoint and a self-adjoint part, 



N = M - S = 



p- x v(VV 2 V) 







(Vu>V) Wc 



(32) 



e(k B T) (%V 2 ) J y [-{VcYV 

where u> is the antisymmetric vorticity tensor, uij^ = dv^/drj — dvj/drk, and we used the vector identity 

(V x v) x v = —v ■ Vv + V 

Even though by skew symmetry the top right sub-block of S is nonzero, there is no coupling of concentration back in 
the velocity equation because 



dH 

dc 



Vc 



dHc 
dc 



Vc = VH C 



is a gradient of a scalar and is eliminated by the projection. The velocity equation therefore remains of the form (30 1 
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B. Spatial Discretization 



For a detailed description of the spatial discretization of ( 30|31 1 that we employ we refer the reader to Ref. [25] , 
The discretization of the velocity equation is based on a staggered or MAC grid [60] in which the component of 
velocity along a given dimension is discretized on a uniform grid that is shifted by half a grid spacing along that 
dimension. Following the stochastic methods of lines that we used for the Burgers equation in Section II B| the spatial 
discretization of ([I]) leads to a system of SODEs of the form 



dv 

~dt 



S v (v) v + vL v v 



2vk B T 
pAV 



D W W V (t) 



(33) 



where S v (v) denotes a discretization of the advective operator — (v ■ V), AV is the volume of a hydrodynamic 
cell, and W„ (i) is a collection of white-noise processes [25]. Here D w a tensor divergence operator that applies 
the conservative discrete vector divergence operator D independently for each coordinate, L v is a discrete (vector) 
Laplacian, and P = I - G (DG)~ D is a discrete projection operator, where G is a discrete scalar gradient operator. 
The imposed periodic, no-slip or free-slip boundary conditions are encoded in the specific forms of the discrete 
difference operators near the boundaries of the domain. 

Let us first focus on creeping Stokes flow, where the advective term v-'Vv is neglected. Following the same procedure 



as we employed for the fluctuating Burgers equation in Section II B the spatial discretization is constructed to obey 
a discrete fluctuation-dissipation balance principle. This relies on several key properties of the staggered difference 
operators. Importantly, the discrete gradient and divergence operators obey the duality relation G = —D*, just as 
the continuum operators. The resulting scalar Laplacian L s = DG = DD* is the standard (2d + l)-point discrete 
Laplacian, which is also applied to each (staggered) component of the velocity to form L v = D w (D w ) (see Ref. 
[2"5] for a discussion of modifications near physical boundaries). Because of the duality between D and G the MAC 
projection is self-adjoint, P* = P, and idempotent, P 2 = P, just like the continuum projection operator. From these 
properties and Eq. (20) it follows (see Appendix in Ref. [35] for details) that the equilibrium covariance of the 
fluctuating velocities is 



(vv ) = 



k B T . 
pAV 



(34) 



This means that when an equilibrium snapshot of the velocity is expressed in any orthonormal basis for the subspace 
of discretely divergence-free vector fields, the coefficients are i.i.d. Gaussian random variables with mean zero and 
variance p~ l ksT / AV . This is the expression of discrete fluctuation-dissipation balance for the case of incompressible 
flow. 

The addition of the nonlinear advective term does not affect the discrete fluctuation-dissipation balance since the 
advective term is skew-adjoint, S* v = —S Vl just like the discretization (12) of the term uu x for the one-dimensional 
case. Specifically, in two dimensions, for a given u such that Du = 0, the spatial discretization of the advective term 
described in Refs. i25j.42! leads to 



[S v (u) v 



(x) 



(4Ax)~ 
(4AyY 



.(*) 



,(*) 



( (x) i (x) \ (x) ( (x) 



This can easily be shown to be a skew-adjoint discretization, [S v (it) v] ■ w = — [S v (u) w) ■ v, for either periodic, 
free-slip and no-slip conditions (or any combination thereof). Furthermore, this discretization leads to Hamiltonian 
dynamics for inviscid flow, i.e., the phase-space flow generated by the advective term is incompressible, 



d_ 

dv 



S v (v) = 0. 



A good temporal integrator should preserve this special structure of the equations and reproduce the correct velocity 
fluctuations for reasonably large time step sizes. 

Note that the addition of a passively-advected scalar field poses no additional difficulties if c is discretized on the 



regular (non-staggered) grid underlying the (staggered) velocity grid. Specifically, the spatial discretization of (31) 
that we employ is a scalar equivalent of (33), 



dc 
dt 



S c (v) c + X L c c + 2 DW C (t) . 



(35) 
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where S c (v) is a cell-centered conservative and skew-adjoint discretization of the advection operator — (v ■ V) 
[S c (u) c] iJ = - (2A X y 1 (u^ .c i+hj - u\ x \ jCi _ ltj j - (2A y y 1 (u ( ^ +1 _c^ +1 - u ( ^^,c id ^ , 



and L c is equivalent to L s = DG except near boundaries. The semi-discrete equation ( 35 ) obeys a discrete fluctuation- 
dissipation balance principle |25j . Specifically, at thermodynamic equilibrium the fluctuations in the concentration 
are Gaussian with covariance 

(36) 



and also uncorrelated with the velocity fluctuations. 

The coupled velocity-concentration equation system of SODEs ( 33|35 ) can be written in the generic Langevin form 
([3]). Note that the concentration-dependent term in the velocity equation that ought to be included to preserve the 
skew-symmetry of the non-diffusive terms can be written in the form of a projected discre te gradient of a scalar, 
¥G (c 2 /2), which vanishes identically. This shows that in two dimensions the system (33 35) is time reversible with 
respect to the equilibrium Gibbs-Boltzmann distribution with the separable discrete Hamiltonian 



k B TAV ^ n 

'•<,.<• 



2c 



1,3 



Note that the Gibbs-Boltzmann distribution is constrained to the linear subspace of discretely divergence-free vector 
fields, 



(Dv) 



Ax' 



Ay- 1 



(y) 
v . 



0. 



and with periodic boundaries both the average momentum and average concentration are conserved by the dynamics. 



C. Temporal Discretization 



In our initial discussion of temporal integration schemes for ( 33 1 we will neglect the advective term and focus on 



creeping Stokes flow, thus avoiding technical details while preserving the essential features of the problem. There is a 
vast literature on deterministic temporal integration of the incompressible Navier-Stokes equations, and, in particular, 
the handling of the V-/T term. One of the most popular class of methods are splitting or projection methods, such as 
the prototype projected Eulcr-Maruyama method, 



= v 



vAt¥L v v n + (2vAt)* FD W WZ, 



(37) 



where we set p~ 1 ksT/ AV = 1 for simplicity and W™ are i.i.d. standard normal random variates generated indepen- 
dently at each time step. Note that in practice, due to roundoff errors and the use of inexact Poisson solvers in the 
projection operation, it is preferable to apply P to v" as well. 

As explained in Appendix [C] (see also Appendix B in Ref. |2S]), the iteration (37) gives a steady-state covariance 
that is a first-order accurate approximation to the continuum result ( 34 ) , 



C v = (v n (v n )*) = V + At AC V + O (At 2 ) . 

In Appendix [C] w e consider approximate projection methods [61U62] and find that they do not satisfy this requirement. 
The scheme (37) can be seen as a direct application of the Euler-Maruyama method to (33). Note that any purely 



explicit scheme, including the RK3 scheme described in Section IIIC can be applied to (33) by simply performing a 
projection operation after every stage of the scheme. 

Semi-implicit schemes can also be applied to (33). As a prototype example, let us consider the implicit midpoint 



method (22 1 



v n+l = v r, 



;At 



■ FL V (v n + v 



{2vAty¥D w Wl. 



At first sight, it appears that solving (38) requires the application of [J — (vAt/2)VL v ] 
see that solving (38) is equivalent to solving the following linear system for the velocity v r 



(38) 

However, it is not hard to 
fl and the pressure 7r™ + 3 , 



vAt 



,n+l 



At Gtt" + 3 
Dv n+1 



vAt 



L v v n + {2vAt) - D w W n v 



(39) 
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This coupled velocity-pressure Stokes linear system can be solved efficiently even in the presence of non-periodic 
boundaries by using a preconditioned Krylov iterative solver, as described in detail in Ref. |63j . The scheme ( 38 ) 
reproduces the static covariance of the velocity fluctuations exactly, C„ = P, for any time step. A more intuitive 



approach to analysing the scheme ( 39 ) based on computing the modes of the spatial discretization is described in 
Appendix ID] 



Having illustrated how the implicit midpoint rule (22) can be applied to the time-dependent Stokes equations, it 



is simple to modify the implicit-explicit predictor-corrector schemes described in Section |III D| to account for incom- 
pressibility. In Ref. [25] the implicit trapezoidal predictor-corrector scheme (28) was used to solve ( 30|31 1, but only 
test ed in an essentially linearized context. Here we also consider the implicit midpoint predictor-corrector scheme 
(29), for which the predictor stage consists of solving the linear system for and c 2 , 
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and similarly for the corrector stage. The concentration equation is decoupled from the velocity equation and can 
be solved using standard techniques, e.g., multigrid methods. Note that in practice it is better to rewrite the linear 
systems in terms of the increments and c n+2 — c™. This is because for weak fluctuations the terms involving 

the identity matrix may dominate the right-hand side and compromise the accuracy of the linear solvers, unless special 
care is taken in choosing the termination criteria for the iterative linear solvers. 



SIMULATIONS 



In this section we numerically study the behavior of the implicit-explicit trapezoidal ( 28 1 and midpoint ( 29 1 schemes 
on the fluctuating Burgers and Navier-Stokes equations. We focus here on the behavior of the equilibrium distribution 
of the fluctuating fields for large fluctuations and large time step sizes. The discrete spectrum of the equilibrium 
fluctuations is one of the most important properties of a numerical scheme for long-time simulations. 

Our spatial discretizations were constructed to obey a discrete fluctuation-dissipation balance principle, which 
means that for sufficiently small time steps the numerical schemes will produce the correct equilibrium fluctuations 



even when the nonlinear terms are important. In the absence of advection, the schemes (28) and (29) reduce to the 
implicit midpoint rule, which was designed to produce the correct equilibrium fluctuations for any time step. It is 
not a priori obvious how increasing the time step size affects the long-time behavior of schemes in the presence of the 
nonlinear advective terms, which are treated explicitly. We study this question numerically in this section. 

In fluctuating hydrodynamics, the magnitude of the equilibrium fluctuations is controlled by the degree of coarse 
graining, more specifically, by the average number of particles N p (microscopic degrees of freedom) per hydrodynamic 
cell (macroscopic degree of freedom) . In particular, the law of large numbers suggests that the Gaussian fluctuations 
at equilibrium have a variance inversely proportional to N p ~ AV. In order to model the effect of the degree of coarse 
graining we can introduce a parameter e ~ N^ 1 that measures the strength of the fluctuations, with e ~ 1 indicating 
very strong fluctuations, i.e., minimal coarse graining. While we cannot expect Markovian SPDE models to be a good 
approximation to reality in the absence of coarse-graining, from a numerical analysis perspective it is important to 
understand how robust the numerical schemes are to increased magnitude of the fluctuations. We note that the types 
of schemes we use here may have uses in other fields such as turbulence, where reproducing the correct spectrum of 
fluctuations is also important. 



A. Fluctuating Burgers Equation 

We first turn our attention to the fluctuating Burgers equation. For simplicity, we take c = 1 and consider a periodic 
system with zero total momentum (no macroscopic advection), 

d t u + ud x u = v d xx u + (2ev) 2 d x Z. 



When the spatial discretization ( 13 ) is used, the coarse-grained velocities Uj have Gaussian equilibrium fluctuations 
with mean zero and covariance 
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The equilibrium magnitude of the velocities is therefore |u| ~ y^T/Ax, which is a measure of the typical magnitude 
of the advection speed. 

1. Dimensionless Numbers 

In deterministic fluid dynamics, the dimensionless number that describes how well advection is resolved by the time 
step size is the advective CFL (Courant-Friedrichs-Levy) number 

M At A 1 . 3 

a = w Ate^Ax~^. 

Ax 

Traditional wisdom says that for schemes that handle advection explicitly At should be chosen such that a < 1. The 
presence of diffusion, however, stabilizes numerical schemes by adding dissipation and introduces the viscous CFL 
number (3 and the cell Reynolds number r, 

vAt a \u\ Ax i i A i 

/?=— , r = - = J-J «e3i/- x Ax2. 

Ax z p v 

Note that the cell Reynolds number measures the relative importance of the nonlinear term (advection) versus the 
linear term (diffusion), and is independent of At. In the deterministic setting, standard von Neumann stability analysis 
for the Eulcr scheme applied to the advection-diffusion equation suggests that our centered discretization of advection 
is stable if a 2 /2 < /3 < 1/2. This suggests that the dimensionless number 

a 2 ,„ Id 2 At eAt 
p v vl\x 

may be important in controlling the behavior of our spatio-temporal discretizations. 

For weak fluctuations (e <C 1) or strong dissipation (r -C 1), both of which are typically true for realistic fluids at 
small scales and sufficient levels of coarse-graining (N p ^> 1), the accuracy is controlled by the viscous CFL number 
P. In particular, if j3 -C 1 we can be confident that the numerical scheme resolves the dynamics of the fluctuations 
accurately. However, running with small viscous CFL numbers is often impractical. Our semi-implicit schemes are 
designed to be stable and also to correctly reproduce the equilibrium fluctuations even for large time steps, /3 > 1, 
at least as long as a < 1. An important question, which is difficult to analyze with existing analytical tools, is how 
large the time step can be before the explicit handling of the nonlinear advective terms introduces large errors. We 
study this question here by examining the equilibrium fluctuations for strong fluctuations, e > 1, and large time steps, 
P > 1. 

2. Static Structure Factors 

In order to study the behavior of the equilibrium fluctuations we follow the approach used in Ref. |41) . Instead of 
studying the fluctuations in the actual (real space) variables Uj, we study the equilibrium discrete Fourier spectrum 
of the fluctuating variables, defined as 

S K = Ne~ x Ax (ukK) ■ 
Here the discrete Fourier transform is defined as 

N 
3=0 

where < k < [N/2\ is the waveindex and Ak = 2itk/N < it is the dimensionless wavenumber. The quantity S K 
is a dimensionless discrete version of what is called the static structure factor S(k) in the physics literature, where 
k = Ak/ Ax is the physical wavenumber. For fluctuating hydrodynamics, at thermodynamic equilibrium 

S K = 1 for all k ^ 0, 



which is a re-statement of the discrete fluctuation-dissipation balance principle. 
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Figure 1: Static structure factor S K for a periodic system of 256 cells at thermodynamic equilibrium. (Left panel) Midpoint 
scheme (291 with Wi = 1/2 and non-Hamiltonian advection (401 for time step size /3 = 0.4 for weak fluctuations e = 0.1 
(a ~ 0.13), medium fluctuations e = 1 (a = 0.4) and strong fluctuations e = 4 (a = 0.8). For e = 4, reducing the time step 
to p = 0.1 (a = 0.2) does not significantly improve the accuracy and even if At — > the wrong spectrum of fluctuations is 
obtained. Switching to the Hamiltonian discretization (|12| | significantly lowers the errors. (Right panel) Comparison between 
the trapezoidal scheme (281 and the midpoint scheme ( |29[ ) with wi = 1/2 and with Wi = 1/4, for e = 4 and large time step size 
P = 0.5 (q = 1) and small time step f3 = 0.125 (a = 0.25). Advection is discretized in a Hamiltonian manner. 



For periodic systems, due to translational invariance, the quantity S K contains the same statistics about the equi- 
librium fluctuations as the N x N covariance matrix Cj t f = (ujUf). The advantage of using the Fourier description 
is that it illustrates the behavior at different physical length scales. It is expected that any numerical scheme will 
produce some artifacts at the largest wavenumbers because of the strong corrections due to the discretization; however, 
small wavenumbers, Afc <C 1, ought to have much smaller errors because they evolve over time scales and length scales 
much larger than the discretization step sizes Ax and At. A scheme or choice of time step size that produces a discrete 
structure factor S K much different from unity at small wavenumbers must be rejected as unphysical. It is important 
to emphasize, however, that getting a good equilibrium static spectrum for the fluctuations is not a guarantee that a 
scheme accurately models the dynamics of the fluctuations. 



3. Numerical Results 



In Fig. [I] we show numerical results for the equilibrium structure factor S K for a periodic system with v = 1 and 
Air = 1 and zero total momentum, (uj) — 0. To illustrate the importance of using a Hamiltonian discretization of 
advection in the nonlinear setting, we consider a scheme where the advective term uu x is handled using the conservative 
but non-Hamiltonian discretization 





instead of the conservative Hamiltonian discretization jl2| ). We recall that the correct answer is S K = 1 for all 
wavenumbers. The results in the left panel of Fig. [I] illustrate that for weak fluctuations (i.e., nearly linear equa- 
tions), the correct spectrum is obtained. However, for strong fluctuations (i.e., nonlinear equations), e = 4, the 
non-Hamiltonian scheme produces the wrong static spectrum of fluctuations at small wavenumbers, and reducing At 
does not help. 

On the other hand, the Hamiltonian discretization gives small errors in the spectrum even for the larger time step 
size. In the right panel of Fig. [I] we zoom in to show the magnitude and form of the errors in the structure factors 



for the implicit-explicit trapezoidal scheme (28) and for the midpoint scheme (291 with Wi = 1/2 and with w\ = 1/4. 
We see that all three schemes show similarly small errors in the spectrum. Reducing At by a factor of four makes the 
error statistically insignificant. We have verified that the errors are of second order in the time step size At for all 
three schemes. 

In Fig. [2] we show numerical results for the static structure factor in the case of weaker fluctuations, e = 0.1 and 
e = 0.01, but large time step size. This is a typical scenario for fluctuating hydrodynamics in practice, since for 
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Figure 2: Static structure factor S K at thermodynamic equilibrium in the case of weaker fluctuations and large time step sizes, 
for the trapezoidal scheme (281 and for the midpoint scheme (29 1 with Wi = 1/2 and with Wi = 1/4. (Left panel) Moderate 
fluctuations, e = 0.1, for time step size /3 = 10, a m 3.2, 7 = 1. (Right panel) Weak fluctuations, e — 0.01 and time time step 
size, /3 = 100, a = 10, 7 = 1. 



reasonable degree of coarse graining the fluctuations would be small and the behavior of the equations would be close 
to that of the linearized equations. In the absence of the advective nonlinearity our schemes are stable for arbitrary 
viscous CFL number /3. Under equilibrium conditions the semi-implicit predictor-corrector schemes we consider are 
observed to be stable up to rather large At as measured in the advective CFL a, and even the dimensionless number 
7 = a 2 / (3. However, for sufficiently large At the nonlinearities are expected to play some role, and one cannot expect 
to be able to increase the time step size up to the stability limit and still maintain reasonable accuracy. The results in 
Fig. [2] show that for large At there appear significant artifacts in the static structure factor for the trapezoidal scheme 
p8| and for the midpoint scheme (29) with w\ = 1/4. While the magnitude of the errors is small, the problematic 
observation is that the errors have a peak at the smallest wavenumbers, where we expect schemes to most closely 
mimic the continuum equations. 

The above observations lead us to select the midpoint scheme (29) with wi = 1/2 as the most robust temporal 
integrator for the fluctuating Burgers equation. At the same time, we should recognize that the best choice of scheme 
will depend on the quantity of interest and the particular problem under consideration. All schemes are observed to 
produce equilibrium fluctuations that are rather robust under the presence of strong non-linearities if a Hamiltonian 
advection of discretization is employed. In the next section we confirm that this conclusion also holds for the fluctuating 
Navier-Stokes equation. 



B. Fluctuating Navier-Stokes Equation with Passive Tracer 



We now turn our attention to the fluctuating Navier-Stokes equations with a passively-advected scalar, Eqs. ( 30|31 1 



Here we set p~ x ksT = e so that e controls the magnitude of both the velocity and the concentration fluctuations 

The implementation of our spatio-temporal discretizations and their accuracy in the linearized setting (weak fluctu- 
ations) is discussed in more detail in Rcf . 25 . Our implementation is integrated into the the IB AMR software frame- 
work 64J, an open-source library for developing fluid-structure interaction models that use the immersed boundary 
method. Based on our experience with the fluctuating Burgers equation, we focus our attention on the implicit-explicit 



trapezoidal scheme (28) and the midpoint scheme (29) with W\ = 1/2. For simplicity, here we focus on two spatial 
dimensions and periodic boundary conditions, but we wish to emphasize that our formulation, numerical schemes, 
and implementation apply to three spatial dimensions and no-slip or free-slip boundaries as well. 

For simplicity, in our numerical tests we set Ax = Ay = 1. The diffusion coefficients for momentum and concentra- 
tion are set to v = 1 and \ — 0.25, and a grid of size 64 x 64 is employed. The same dimensionless numbers as for the 
Burgers equation apply, with the difference that there is a separate diffusive CFL for the concentration, j3 c — xP/ v i 
and therefore the cell Peclet number r c — \u\ Ax/\ is four times larger than the cell Reynolds number r. 
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1. Static Structure Factors 

The equilibrium fluctuations in velocity and concentration are characterized by the static structure factors, which are 
the equilibrium average of the discrete Fourier spectrum of the fluctuating velocities and concentrations. Concentration 
fluctuations are characterized via 



where AV^ = Ax Ay is the volume of the hydrodynamic cells, N c is the number of hydrodynamic cells, and k, = (k x , k v ) 
is the waveindex. For the velocity fluctuations, we calculate the spectrum of the fluctuations of a variable related to 
vorticity [2"5] . 

s£ n > = N^AV (fi K tit 



where f2 K is obtained from the discrete Fourier spectrum of the velocity components as 

^K, = k {k x Vy kyVx) , 



and k — ^Jk x + ky is the wavenumber. Note that fully characterizes the covariance of the velocity fluctuations 

since our scheme ensures the velocity is discretely divergence free at all times and k^ 1 (—k y , k x ) spans the subspace of 
divergence-free velocities in Fourier space. For staggered variables the shift between the corresponding grids should be 
taken into account as a phase shift in Fourier space, for example, exp (k x Ax/2) for v x . Additionally, the wavenumber 
k = (k x ,k y ) should be replaced by the effective wavenumber k that takes into account the centered discretization of 
the projection operator, for example, 

~ exp (ik x Ax/2) — exp (— ik x Ax/2) sin (k x Ax/2) 
x ~ iAx ~ x (k x Ax/2) ' ' 

One can additionally define and measure the cross-correlation between concentration and velocity fluctuations via the 
cross-correlation static structure factor 

= N^AV (c R fr K 



which is in general a complex number. 

Discrete fluctuation-dissipation balance in the coupled velocity-concentration equations requires that = = 
1 and St v) = for all nonzero wavenumbers. Deviations from these values indicate a violation of discrete fluctuation- 
dissipation balance and can be used to numerically assess the behavior of the schemes in the nonlinear setting, as we 
do next. 

2. Numerical Results 

In Fig. [3] we show numerical results for the spectrum of the fluctuations in the solenoidal modes of velocity 
and concentration for strong fluctuations, e = 4. We see that, just as for the fluctuating Burgers equation, both the 
trapezoidal and the midpoint scheme show artifacts in the spectra, especially for concentration and for the trapezoidal 
scheme. The cross-correlation is found to be small and difficult to measure due to large statistical errors. 

We have verified that as the time step is reduced, both schemes give the correct spectrum even for strong fluctuations 
(i.e., strong nonlinearities). In Fig. 4]we show the average error in the equilibrium spectrum (\S K — 1|) for vorticity 
and concentration. The second-order weak accuracy of the error is clearly seen in Fig. [2] for both velocity and 
concentration and for both the midpoint and the trapezoidal scheme. 

In Fig. [5] we show the spectrum of concentration fluctuations for the case of weak fluctuations, e = 0.01, and 
large time step size, viscous CFL j3 — 50 and diffusive CFL f3 c — 12.5, and advective CFL a = 5. We see a large error 
for small wavenumbers for the trapezoidal scheme. A similar but weaker artifact is seen for the velocity spectrum 
as well, indicating that the trapezoidal scheme violates fluctuation-dissipation balance at small wavenumbers for 
large At. The midpoint scheme is seen to be much more accurate for both Sk^ and . These investigations confirm 



that the midpoint scheme (29) with u>\ = 1/2 is the more robust temporal integrator for fluctuating hydrodynamics. 
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Figure 3: Discrete structure factors S« ' (top) and (bottom) for the midpoint scheme (291 with wi = 1/2 on the left and 
for the trapezoidal scheme ( 28 \ on the right, for the case of large fluctuations e = 4. ( Top row) Velocity spectra for a time step 



size /3 — 0.5, a = 1. {Bottom row) concentration spectra for /3 C = xP/ v — 0.0625 and a = 0.5. 



VI. CONCLUSIONS 



Discrete and continuum Langevin models are often used as coarse-grained models for the behavior of materials 
at mesoscopic scales. These models include the effects of thermal fluctuations via white-noise stochastic forcing 
terms chosen in a way that ensures fluctuation-dissipation balance. This means that at thermodynamic equilibrium 
the dynamics is time reversible (i.e., in detailed balance) with respect to the Gibbs-Boltzmann distribution. For 
continuum models this is usually only true formally and a more precise interpretation of the equations requires 
introducing a spatial discretization and truncating the continuum models at a scale well-separated from the molecular 
scale. Here we focused on fluctuating hydrodynamics, specifically, we considered numerical methods for solving the 
fluctuating Burgers equation in one dimension and the fluctuating Navier-Stokes equations in two dimensions. In these 
equations, fluctuation-dissipation balance is obtained from the balance of the dissipative (self- adjoint) diffusive terms 
and the stochastic forcing. The advection terms are non-dissipative (skew-adjoint) and do not affect the equilibrium 
distribution. 

Spatio-temporal discretizations do not necessarily preserve the properties of the continuum equations, notably, they 
may not obey a discrete fluctuation-dissipation principle. This means that their equilibrium distribution is not a 
discrete form of the Gibbs-Boltzmann distribution. This sort of unphysical behavior can be avoided by carefully 
constructing the spatial discretization to obey a discrete-fluctuation principle with respect to a target discrete Gibbs- 
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Figure 4: Average error in the static spectra for several time step sizes for the trapezoidal scheme 0288 and for the midpoint 
scheme (29 1 with wi = 1/2, for the case of large fluctuations e = 4. First- and second-order error trends are indicated, showing 

the second-order asymptotic accuracy. Error bars are comparable to the symbol size and not shown. (Left panel) 
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Figure 5: Discrete structure factor for concentration for the midpoint scheme (29 1 with Wi = 1/2 on the left and for the 



trapezoidal scheme (281 on the right, for the case of weak fluctuations, e = 0.01, and large time step size /3 = 50, /3 C = 12.5, 
a = 5. 



Boltzmann distribution. In this way, a coarse-grained semi-discrete model is obtained that obeys the principles of 
statistical mechanics, namely, obeys detailed balance at thermodynamic equilibrium. We showed explicitly how one 
can construct such spatial discretizations by first understanding the properties of the continuum model (even if only at 
a formal level) and then maintaining those properties in the spatial discretization. For fluctuating hydrodynamics, this 
means that certain relations between the discrete divergence, gradient and Laplacian operators must be preserved, such 
as the fact that the divergence and gradient operators are negative adjoints of each other. We showed how to construct 
spatial discretizations that obey a discrete fluctuation-dissipation principle for both the fluctuating Burgers and the 
fluctuating Navier-Stokes with the addition of a passively-advected scalar. These results were mostly a summary of 
discretizations previously constructed in somewhat disjoint bodies of literature, and are a stochastic equivalent of a 
method-of-lines approach for deterministic fluid dynamics [24] , 

The novel challenge that we tackled here is the construction of temporal integrators that are weakly second-order 
accurate for additive noise, and preserve the fluctuation-dissipation balance reasonably accurately. It is possible to 
construct temporal integrators that strictly preserve the Gibbs-Boltzmann distribution at equilibrium by combining 
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a Metropolis-Hastings rejection procedure with a more classical temporal integrator. However, effective use of such 
Metropolization [53] requires first constructing a classical weak integrator that produces a good approximation to 
the equilibrium distribution for time step sizes that are only limited by the deterministic dynamics. For separable 
Langevin equations, it is well-known how to construct symplectic integrators that preserve dynamical invariants, 
including the equilibrium distribution, very robustly [23] . One of the simplest second-order methods of this type 
suitable for fluctuating hydrodynamics is the implicit midpoint rule. However, an implicit discretization of the 
nonlinear advective terms is arguably impractical in three dimensions. Many existing deterministic schemes in fluid 
dynamics use a mixed implicit-explicit approach in which the diffusive terms are handled implicitly and the advective 
terms are handled explicitly. We showed how to construct predictor-corrector second-order Runge-Kutta schemes 
based on the implicit midpoint rule (Crank-Nicolson method) for the diffusive terms, and either an explicit midpoint 
or an explicit trapezoidal rule for the advective terms. These schemes have the remarkable property that in the 
linearized setting (weak fluctuations) they exactly preserve the Gaussian approximation to the Gibbs-Boltzmann 
distribution. We obtained the conditions for second-order weak accuracy of a rather general class of two-stage Runge- 
Kutta methods, and proposed several specific schemes. Through numerical investigations, we found the scheme (29 1, 
which combines the implicit and explicit midpoint rules, was very robust in reproducing the correct spectrum of 
fluctuations, even for strong fluctuations (strong nonlinearities) and for large time step sizes. 

A remaining challenge is how to control not just the static distribution of the equilibrium fluctuations but also 
to control the accuracy of the dynamics of the fluctuating fields. This is particularly challenging for the case in 
which there are several physical processes with large separation of time scales, such as, for example, the diffusion 
of momentum and of mass. Notably, for many realistic fluids the Schmidt number S c — v/x is on the order of 
10 2 — 10 3 , making the sort of integrators we discussed impractical because the time step would be severely limited 
by the viscosity v and not by the diffusion coefficient x- Temporal integrators for this sort of multiscale, or more 
appropriately, manyscale system of SPDEs, will be the subject of future investigations. Future work should also 
explore Metropolization |53| strategies for predictor-corrector schemes for fluctuating hydrodynamics, particularly for 
the case of multiplicative noise. 

We focused here on constructing temporal integrators with weak accuracy of order two for additive equations. This 
choice is guided in part by physical intuition that has yet to be justified more rigorously. Specifically, fluctuating 
hydrodynamics is typically applied at length and time scales where fluctuations are weak and the dynamics is in a 
nearly linearized regime. Notably, in many cases of interest the noise is essentially additive or can be well-approximated 
as additive. Weak accuracy is emphasized because it is more relevant in the types of Monte Carlo applications we are 
interested in. While the classical definition of strong order of accuracy is likely too strong for Monte Carlo simulation, 
some notion of pathwise convergence is important in reproducing the statistics of typical paths, such as, for example, 
rare transition statistics. Recently, error metrics other than weak and strong error have also been considered. Notably, 
Ref. [55] defines the numerical error as the difference between the joint PDF of the computed values at all time steps 
and the joint PDF of the exact process at the same times, and constructs Runge-Kutta integrators accurate in that 
metric. Such approaches may be fruitfully applied in fluctuating hydrodynamics in the future, especially in the context 
of modeling rare events. 
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Appendix A: Weak Semi-Implicit Predictor-Corrector Schemes 



We consider implicit-explicit stochastic Runge-Kutta schemes to solve the system of Ito SODEs 

dx = L(x)xdt + g(x)dt + ad x ■ [K (x)K*(x)} dt + K{x) dS{t), 

where L(x) is a linear operator that represents the implicitly-treated part of the dynamics. Here B is a collection of 
independent Brownian motions (Wiener processes), with W = dB/dt denoting a collection of white noise processes. 
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For increased generality we allow for multiplicative noise and include a divergence term proportional to some constant 
a, as may be required to ensure fluctuation-dissipation balance in the case of multiplicative noise. For example, 
a = 1/2 gives the spurious or thermal drift term in the generic Langevin equation |3]) if the non-dissipative part of 
the dynamics is Hamiltonian. We develop a method that is weakly second order accurate for additive noise and first 
order accurate for multiplicative noise, requires solving only linear systems involving the matrix L n = L(x n ), and 
does not require evaluating the divergence of the mobility matrix M = KK* . 

We employ a two-stage (predictor-corrector) Runge-Kutta method, where at time step n the predictor estimates a 
first-order accurate solution at an intermediate time nAt + w 2 At, 



x = x n + (w 2 - W!)AtL n x n + Wl AtL n x + w 2 Atg n + (w 2 At)* JTWJ, 

and the corrector evaluates the solution at time (n + I) At, 

x n+i = x n + (i- W3 - W4 )AtL n x n + w 3 AtL n x + w 4 AtL n x n+1 
+ w 5 At(L - L n )x + w 5 Atg + (1 - w b )Atg n 
+ Ai2 [(1 - w 6 )I + w 6 M (M")- 1 } K n \wjW'l + {1 - w 2 y Wl 



(Al) 



(A2) 



The handling of the multiplicative noise term is inspired by the so-called kinetic interpretation of the stochastic 
integral [36] and the well-known Fixman method for Brownian dynamics |66| . In the above discretization, the standard 

normal variates W™ correspond to the increment of the underlying Wiener processes over the time interval w 2 At, 

i 

3 (nAt + w 2 At) — 3 (nAt) = (w 2 At) 2 W™, while the normal variates W 1 ^ correspond to the independent increment 
over the remainder of the time step, 3 ((n + 1) At) - 3 (nAt + w 2 At) = ((1 - w 2 ) A<) 3 W 2 . 



1. Additive Noise 



For additive noise, K(x) = K, we would like to achieve second-order weak accuracy. A second-order integrator that 
uses derivatives is provided by the weak Taylor series (161, in indicial notation with the implied summation notation, 



x a (t + At) 



x a (t) 
At 2 



AtL a pXp 
At 2 



Atg n a + K a pAB 



At 



At 2 



K ie K Se (x^dsL^ + 29 7 i2 6 + d^dsgl) + O (At*) , 



(A3) 



1/2 

where A3 = 3 ((n + 1) At) — 3 (nAt) is the Wiener increment, in law, A3 — (At) 1 W where W is a collection of 
i.i.d. standard normal random variables. 

To prove second-order order weak accuracy for the derivative-free Runge-Kutta scheme, we need to match the 
first 5 moments of the numerical increment A™ = — iJJ to the moments of the true in crement A a (nAt) = 

x a ((n + 1) At) — x a (nAt), up to order At 2 [51]. The increment of the RK method ( Al|A2 ) is (after recursively 
substituting in for x n+1 , and Taylor expanding certain terms) 



„n+l 



= x n a + AtL^pXp + Atgl + At*K afj [wi {W?) p + (1 - w 2 y (W?) p 
(At 2 L^ + At 2 g^) [w 2 w 5 x n pd 7 L 7 >p + (w 2 w 3 + w 4 )L^ + w 2 w 5 d ig ™ 
At 3 / 2 K. 



■ye 



w 5 wixp (d^) (WT) e +Uw 3 + w 4 )wi (W?) e +w 4 (l- w 2 y (Wf), 



0:7 



■W5«>Z&g2)(W?) t 



W 2 W 5 



At 2 K 7t K Se (xpd^dsL^p + + d 7 d s g2) + O Ut 5 / 2 



Comparing the first moments of the increments, E (A™) — E (A a (nAt)) — O (At 3 ) if w 2 w$ = 1 and w 2 ws + w 4 = h. 
The difference in the second moments is also O (Ai 3 ) under the same conditions, as is the difference in the third and 
fourth moments. The fifth moments are already of O (At 3 ) . 
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2. Multiplicative Noise 

For multiplicative noise, we only aim to achieve first-order weak accuracy. For this, we need to match the first three 
moments of the increment up to first order in At [51]. With a variable K(x), the difference in the first moments 
obtains additional terms 

(a - w 2 w 6 )Atd J (K aP K^) + 0{At 2 ) 

Therefore, if we have w 2 wq = a, then our method will be consistent with the Ito SDE and reproduce the correct 
thermal drift without explicitly evaluating d x ■ M(x). Note however that this type o f "Fi xman" approach requires 



solving one more linear system per time step because of the appearance of (M n ) 1 in (A2). 

Despite recent progress in mathematical understanding [67] , multiplicative noise has unclear physical relevance. 
Fluctuating hydrodynamics is a coarse-grained model in which the variables are averages over many microscopic 
degrees of freedom and therefore the fluctuations are weak @HJ EH] . To leading order in the magnitude of the 
fluctuations, the equations are linear and the noise can be considered additive (though potentially time-dependent). 
We speculate that multiplicative noise is a very weak perturbation and its effects on the larger scales of the flow, if 
any, can be sufficiently accurately captured by low-order temporal integrators such as ( A1|A2 ) with the condition 
w 2 wq = a. 



Appendix B: L-Stable Scheme 



In typical fluctuating hydrodynamics applications, for explicit schemes the time step is severely limited not by 
advection but by momentum or heat diffusion, notably, by viscous dissipation. For purely dissipative linear equations, 



implicit handling of momentum diffusion can yield A-stable schemes such as the implicit midpoint scheme (22 1. This 



allows the use of much larger time step size At, at least in principle. If one is interested in steady-state fluctuations, 
the implicit midpoint scheme ( 22 ) gives the correct spectrum of fluctuations for any At (see the Appendix in Ref. [25] 
for a discussion of how to choose a suitable At) . 

However, for time dependent linear problems, only an exponential integrator can reproduce the correct dynamics 
for all modes (wavenumbers) for all time step sizes. The implicit midpoint rule provides a notably bad approximation 
to the exponential decay of correlations for large At, since the Pade (1,1) rational approximation to the exponential 
p3| , exp(— x) f=3 (1 — x/2)/(l + x/2) tends to -1 for x ^> 1 instead of decaying to zero. This leads to oscillatory 
dynamics for the modes that are under-resolved by the large time step size, i.e., for the thermal fluctuations at large 
wavenumbers. A much better approximation to exp(— a;) is provided by rational approximations that decay to zero 
as x — > 00. In numerical analysis jargon this means handling the diffusive fluxes using an L-stable numerical method. 

Let us consider the choice of weights in the general scheme ( 25 ) that yield a scheme that is weakly second-order 
accurate and L-stable in the implicit part of the dynamics. From the conditions for second-order accuracy (26) we 
obtain 



w 3 = 



and from the condition of L stability we obtain 



W4 



w 2 



w 5 



1 

2w 2 



w 1 = 



Wi 
Wi 



which gives the following rational approximation to the exponential decay of the dynamics, 

(l - 2 Wi + 2 u>!) x - 2 (1 - 104) 



exp (—x) 



(w 4 x + 1) [(2 Wi - 1) x - 2 (1 - wi)] ' 



(Bl) 



A reasonable choice of w 4 can be taken to be the one that minimizes the mismatch between the coefficient in front of 
x 3 in the Taylor series expansion of the left and right hand sides [SS], giving w 4 = 1 ± y/2/2. In Fig. [HJwe compare 
the two rational approximations by choosing the plus or minus sign. In the deterministic literature the choice of the 
minus sign has been favored }50l 169] . however, we recommend the plus sign, 



2G 




Crank-Nicolson 
L-stable (plus sign) 
L-stable (minus sign) 



Figure 6: Comparison of exp(— x) with three rational a ppro ximations. The approximation in the Crank-Nicolson scheme (231 
does not decay to zero for large x. The approximation (Bl I decays to zero for u>4 = 1 ± \/2/2, however, only the positive sign 
gives a strictly positive approximation. 



because this gives a strictly positive approximation to the exponential decay of correlations instead of oscillatory 
behavior for the under-resolved modes (large x). 

We are still left with the choice of W2, where the two common choices for W2 would be a mid-point, W2 = 1/2, or an 
end-point, W2 
explicitly (L 



- 1, predictor stage. 
0), the choice W2 - 



From the discussion in Section |IIID 1| we know that when all terms are handled 

It 



1/2 gives third-order accuracy for the static covariance for linear problems 
can also be shown that this choice leads to third-order accuracy of the static covariances in the linearized setting if 
all terms are discretized implicitly (g(x) = 0). This suggests that a better, even if not unique, choice, is to take 
W2 — 1/2, giving our preferred choice of weights for an L-stable predictor-corrector scheme, 



1 



V2 



1 



U> 2 



, w 3 = -(1 + V2), 



1. 



(B2) 



In the linearized setting, this L-stable scheme gives second-order accurate covariances for small time step sizes; however, 
it does not produce the correct spectrum for the fluctuations for large time step sizes, unlike the implicit midpoint 
scheme (22 1. In particular, for large At the L-stable scheme strongly damps the magnitude of the fluctuations of the 



fast (small wavelength or large wavenumber) modes. Therefore, if static covariances are the quantity of interest, the 
implicit midpoint rule should be used instead. 



Appendix C: Approximate Projection Methods 



Here we consider a generalization of the projected Euler-Maruyama scheme ( 37 1 , 

n + vAt L v v n + (2^Ai)5 D w Wl \ , 



v n+l = 



(CI) 



where P an approximation to the discrete projection P, for example, P = I — GL~ X D. Here L p is a discrete pressure 
Laplacian operator that may, in general, be different from L s = DG. For example, with spatial discretizations of 
the incompressible (Navier-) Stokes equations that use cell-centered velocities, L s possesses a non-trivial nullspace 
and the corresponding exact projection methods (in which P — P) suffer from the so-called checkerboard instability. 
Approximate projection methods have been developed to overcome these difficulties of exact cell-centered projection 
methods [61j . One of the simplest approximate projection methods is (CI) with L p being the standard second-order 
Laplacian stencil [55]. For the staggered-grid spatial discretization we employ, however, it is straightforward to invert 
L s and approximate projection methods arc not used in practice. 

The steady-state covariance of the iteration (CI I should be a consistent approximation to the continuum result 



(34). Specifically, we ask that to leading order in the time step size 



-1 ( V «H 



(v n (v n )*) = ¥ + AtAC v + O (At 2 
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Substituting (CI I in this condition and equating the leading-order terms we obtain the condition 



This condition is satisfied for exact projection methods, P = P, but not for approximate projection methods, P ^ 



Assuming the initial condition is discretely divergence-free, for exact projection (CI) is equivalent to (37) 



Appendix D: Mode Analysis 



It is instructive to describe a framework for analyzing schemes such as (39), following the mode analysis used 
to study splitting errors in projection methods in the deterministic context 70| I71j. This analysis can in principle 
produce explicit expressions for the spectrum of velocity fluctuations for the types of schemes we consider here. It 
also illustrates clearly the role of the pressure and, in particular, the difficulties with applying semi- implicit projection 
(splitting) methods in the context of the fluctuating Navier-Stokes equations. 

A mode of the spatially-discretized unforced time-dependent (creeping) Stokes flow equation 

d t v + Gtv = vL v v, s.t. Dv = Q (Dl) 

is an exponentially-decaying solution of the form 

v (t) = v a er at and tt (t) = VTT er at . 

Here a > is the decay rate associated with the spatial mode {vq, ttq}, which is a normalized solution to the 
eigen-problem 

(L v + v~ x at) v Q + Gn Q = and Dv = 0. (D2) 

These modes diagonalize the creeping Stokes flow dynamics and form a complete orthonormal basis for the space of 
divergence-free velocity fields. This can be seen by eliminating the pressure to obtain the classical eigenvalue problem 
in the subspace of discretely divergence-free velocity fields 



L v — L^G (DL~ X G) X DL- x 



v a = -v 1 av . 



In the presence of a stochastic forcing, we can express any solution in a basis formed by the modes {vq, Vq, . . . }, 

»(*)=2«k(*)«oi 

k 

l 

where the mode amplitudes Vk(t) are scalar stochastic processes. The stochastic forcing (2v) 2 D w y\? v {t) in the 
momentum equation can be projected onto Vq to obtain the amplitude of the stochastic forcing for mode k, 

w k (t) = (2v)* («J)*(D„W, (t)), 
which is a scalar white-noise process with covariance 

(w k (t) wl (f )) = lv (v*Y [D w <W„ (t) W: (*')> D* w ] = -lv («*)* L v (v k ) 6 (t - t') , 

where we made use of the discrete fluct uation-dissipation balance between the viscous dissipation and the stochastic 
forcing, L v = D w (D w )* . From (D2) we can express 

-2v («§)* L v = -2v («§)* {Gtvq - v^o-kvf) = 2v (Dv$)* n + 2a k \\v%\\ 2 = 2(7*, 

where we again made use of the duality relation G = —D*. This simple calculation shows that in the mode represen- 
tation the linearized fluctuating Navier-Stokes equation becomes a collection of decoupled scalar Langevin equations 
driven by standard Wiener processes, 

C ^ = -^* + (2 ( 7 fe )'W„(t). (D3) 
at 
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The fluctuation-dissipation balance is most clearly revealed in this representation. 

Temporal discretizations can be analyzed by projecting the numerical solution onto a set of discrete modes. For 
the implicit midpoint discretization ( 39 1 , the modes are decaying solutions of the form 



and 



n+i 



where &k ps o> is the numerical decay rate. The spatial structure of the mode {v%, 7rJJ} is the solution to the discrete 
eigen-problem 



L v + v 



e -a k At _|_ I 



At 



Vi. 



1 



Ur7Tt. 



Comparison to (|D2j) shows that the spatial modes are the same as for the semi-continuum (Dl), and the temporal 
decay rate is second-order accurate in the time step, 



2 

g — CTfc At 



1 



1 -e 
Af 



-CTfcAt 



= ff* 1 - 



12 



O (At 3 



When the stochastic forcing is included, the discrete velocity can be represented in the basis formed by the discrete 
modes just as we did above for the time-continuous equations. In the mode representation the scheme (39) is seen to 
be nothing more than the implicit midpoint method (22 1 applied to the system of decoupled SDEs (D3|. 



The mode analysis reveals that semi-implicit projection (splitting) methods have a significant shortcoming not seen 
for explicit methods. A Crank-Nicolson projection method for (Dl I consists of first solving the following linear system 
for the velocity v n+1 with a time-lagged pressure [72] . 



At Gi 



vAt 



and then projecting the intermediate velocity v n+1 
a linear system for the pressure correction 



to enforce the divergence-free constraint, v 



n+l 



by solving 



v n+i = -n+i _ AtGAvr", s.t. Dv n+1 = 0. 



Repeating the discrete mode calculation reveals that the spatial modes for the above temporal discretization are not 
the same as for the semi-continuum (Dl), specifically, the gradient of pressure term in (D2| is modified by a term 
involving the Laplacian L v . For periodic systems the discrete gradient and vector Laplacian commute, L V G = GL S , 
and modes have the correct spatial structure. However, for non-periodic systems the splitting of the pressure and 
velocity equations introduces a commutator error that leads to the appearance of "spurious" or "parasitic" modes [7T] . 
For deterministic solutions and moderate time step sizes, spatio-temporal smoothness of the solution usually makes 
these commutator errors acceptably small. In the stochastic context, however, all modes are stochastically forced and 
have a non-negligible amplitude, including the parasitic modes. For this reason, we chose to use (39) and solve a 



coupled Stokes linear system for both pressure and velocity, and only use the projection method as a preconditioner 
for the required Krylov solver |63] . We emphasize again that for purely explicit time stepping scheme the spatial 
structure of the modes is preserved and projection methods can be used in the stochastic setting as well. 
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